In this script, I fit simplified models to density data so that I can predict densities on the condition data and pred grid
rm(list = ls())
library(tidyverse); theme_set(theme_light(base_size = 12))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/clean_data/cod_fle_density_models_as_covars_cache/html")
world <- ne_countries(scale = "medium", returnclass = "sf")
# Specify map ranges
ymin = 54; ymax = 58; xmin = 12; xmax = 22
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
#ggplot(swe_coast_proj) + geom_sf()
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")
# Filter to match the data I want to predict on
d <- d %>% filter(year > 1992 & year < 2020 & quarter == 4)
# Calculate standardized variables
d <- d %>%
mutate(depth_sc = (depth - mean(depth))/sd(depth)) %>%
mutate(year = as.integer(year)) %>%
drop_na(depth) %>%
rename("density_cod" = "density") # to fit better with how flounder is named
# Explore data a bit
ggplot(d) +
geom_point(aes(year, density_cod, color = factor(sub_div))) +
facet_wrap(~sub_div)
# Explore data a bit
ggplot(d) +
geom_point(aes(year, density_fle, color = factor(sub_div))) +
facet_wrap(~sub_div)
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
#> Rows: 129346 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
#> Rows: 120107 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- bind_rows(pred_grid1, pred_grid)
# Standardize data with respect to data
pred_grid <- pred_grid %>%
mutate(year = as.integer(year)) %>%
filter(year %in% c(unique(d$year))) %>%
mutate(depth_sc = (depth - mean(d$depth))/sd(d$depth),
temp_sc = (temp - mean(d$temp))/sd(d$temp),
oxy_sc = (oxy - mean(d$oxy))/sd(d$oxy)) %>% # Need to scale these to the mean and sd in the data!
drop_na(oxy, depth, temp)
#> mutate: converted 'year' from double to integer (0 new NA)
#> filter: no rows removed
#> mutate: new variable 'depth_sc' (double) with 6,238 unique values and 0% NA
#> new variable 'temp_sc' (double) with 249,453 unique values and 0% NA
#> new variable 'oxy_sc' (double) with 249,453 unique values and 0% NA
#> drop_na: no rows removed
# Make mesh:
spde <- make_mesh(d, xy_cols = c("X", "Y"),
n_knots = 200,
type = "kmeans", seed = 42)
plot(spde)
# Depth spline + oxy spline
# Cod model
mcod <- sdmTMB(density_cod ~ 0 + as.factor(year) + s(depth_sc),
data = d,
mesh = spde, family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_steps = 1))
tidy(mcod, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 as.factor(year)1993 4.837552 0.4688664 3.918590 5.756513
#> 2 as.factor(year)1994 5.197526 0.4517869 4.312040 6.083012
#> 3 as.factor(year)1995 5.807254 0.4443790 4.936288 6.678221
#> 4 as.factor(year)1996 5.083857 0.4449677 4.211736 5.955977
#> 5 as.factor(year)1997 4.316940 0.4287832 3.476540 5.157339
#> 6 as.factor(year)1998 4.502834 0.4325650 3.655022 5.350645
#> 7 as.factor(year)1999 4.347810 0.4224004 3.519921 5.175700
#> 8 as.factor(year)2000 4.694878 0.4021157 3.906746 5.483011
#> 9 as.factor(year)2001 5.099759 0.3908475 4.333712 5.865806
#> 10 as.factor(year)2002 5.712133 0.3823567 4.962728 6.461539
#> 11 as.factor(year)2003 4.502260 0.4054629 3.707568 5.296953
#> 12 as.factor(year)2004 4.959874 0.4083146 4.159592 5.760156
#> 13 as.factor(year)2005 5.346098 0.3743546 4.612376 6.079819
#> 14 as.factor(year)2006 5.041373 0.3600782 4.335632 5.747113
#> 15 as.factor(year)2007 5.280502 0.3580791 4.578680 5.982324
#> 16 as.factor(year)2008 5.495251 0.3536101 4.802188 6.188314
#> 17 as.factor(year)2009 5.438676 0.3608035 4.731514 6.145838
#> 18 as.factor(year)2010 5.418556 0.3576511 4.717572 6.119539
#> 19 as.factor(year)2011 4.679926 0.3587164 3.976855 5.382997
#> 20 as.factor(year)2012 4.422424 0.3662989 3.704492 5.140357
#> 21 as.factor(year)2013 4.735018 0.3622512 4.025019 5.445017
#> 22 as.factor(year)2014 4.582177 0.3638320 3.869079 5.295275
#> 23 as.factor(year)2015 4.681539 0.3634106 3.969268 5.393811
#> 24 as.factor(year)2016 4.124882 0.3637063 3.412030 4.837733
#> 25 as.factor(year)2017 4.560351 0.3605213 3.853742 5.266960
#> 26 as.factor(year)2018 3.664196 0.3665505 2.945770 4.382621
#> 27 as.factor(year)2019 3.762937 0.4148706 2.949806 4.576069
sanity(mcod)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_tau_O` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_tau_E` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_kappa` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_smooth_sigma` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcod_1 <- run_extra_optimization(mcod, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_1)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> ✓ No gradients with respect to fixed effects are >= 0.001
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcmc_res <- residuals(mcod, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#> Constructing atomic tweedie_logW
#> Constructing atomic invpd
#>
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.006393 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 63.93 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
#> Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
#> Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
#> Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
#> Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 1041.76 seconds (Warm-up)
#> Chain 1: 6.70257 seconds (Sampling)
#> Chain 1: 1048.47 seconds (Total)
#> Chain 1:
qqnorm(mcmc_res);qqline(mcmc_res)
# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mcod_1, "output/mcod.rds")
# Fit flounder model
mfle <- sdmTMB(density_fle ~ 0 + as.factor(year) + s(depth_sc),
data = d,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_steps = 1))
tidy(mfle, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 as.factor(year)1993 3.556486 0.6490165 2.284437 4.828535
#> 2 as.factor(year)1994 3.083573 0.6494050 1.810763 4.356384
#> 3 as.factor(year)1995 4.232611 0.6346263 2.988767 5.476456
#> 4 as.factor(year)1996 3.618856 0.6370869 2.370188 4.867523
#> 5 as.factor(year)1997 3.397782 0.6169767 2.188530 4.607034
#> 6 as.factor(year)1998 2.783891 0.6189851 1.570703 3.997080
#> 7 as.factor(year)1999 2.350747 0.6129248 1.149437 3.552058
#> 8 as.factor(year)2000 3.154737 0.5957873 1.987016 4.322459
#> 9 as.factor(year)2001 3.575064 0.5883076 2.422002 4.728125
#> 10 as.factor(year)2002 4.537331 0.5806431 3.399291 5.675371
#> 11 as.factor(year)2003 3.375738 0.5913574 2.216699 4.534777
#> 12 as.factor(year)2004 3.844600 0.5955993 2.677247 5.011953
#> 13 as.factor(year)2005 3.704495 0.5725463 2.582324 4.826665
#> 14 as.factor(year)2006 3.668069 0.5594159 2.571634 4.764504
#> 15 as.factor(year)2007 3.942359 0.5580117 2.848676 5.036041
#> 16 as.factor(year)2008 4.099267 0.5556642 3.010185 5.188349
#> 17 as.factor(year)2009 4.186882 0.5603307 3.088654 5.285110
#> 18 as.factor(year)2010 4.225439 0.5552249 3.137218 5.313659
#> 19 as.factor(year)2011 3.853699 0.5550766 2.765769 4.941629
#> 20 as.factor(year)2012 3.604305 0.5573405 2.511938 4.696673
#> 21 as.factor(year)2013 3.837236 0.5587290 2.742147 4.932324
#> 22 as.factor(year)2014 3.606003 0.5583835 2.511592 4.700415
#> 23 as.factor(year)2015 3.746447 0.5577288 2.653319 4.839575
#> 24 as.factor(year)2016 3.768052 0.5565157 2.677301 4.858803
#> 25 as.factor(year)2017 3.901198 0.5575261 2.808466 4.993929
#> 26 as.factor(year)2018 3.636791 0.5605948 2.538045 4.735536
#> 27 as.factor(year)2019 3.920627 0.5897827 2.764674 5.076580
sanity(mfle)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_phi` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mfle_1 <- run_extra_optimization(mfle, nlminb_loops = 1, newton_loops = 1)
sanity(mfle_1)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> ✓ No gradients with respect to fixed effects are >= 0.001
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcmc_res_fle <- residuals(mfle_1, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#>
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.006621 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 66.21 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
#> Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
#> Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
#> Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
#> Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 1030.14 seconds (Warm-up)
#> Chain 1: 5.2098 seconds (Sampling)
#> Chain 1: 1035.35 seconds (Total)
#> Chain 1:
#> Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
qqnorm(mcmc_res_fle);qqline(mcmc_res_fle)
# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mfle_1, "output/mfle.rds")
pred_gid (make_pred_grid_utm.R)### Cod
# SD 24
preds_cod24_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 211,032 rows (85%), 38,421 rows remaining
# SD 25
preds_cod25_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 178,983 rows (72%), 70,470 rows remaining
# SD 26
preds_cod26_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 184,410 rows (74%), 65,043 rows remaining
# SD 27
preds_cod27_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 224,505 rows (90%), 24,948 rows remaining
# SD 28
preds_cod28_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 198,882 rows (80%), 50,571 rows remaining
# Now calculate the CPUE index (average)
index24_cod_sim <- get_index_sims(preds_cod24_sim, area = rep(2*2, nrow(preds_cod24_sim)))
index25_cod_sim <- get_index_sims(preds_cod25_sim, area = rep(2*2, nrow(preds_cod25_sim)))
index26_cod_sim <- get_index_sims(preds_cod26_sim, area = rep(2*2, nrow(preds_cod26_sim)))
index27_cod_sim <- get_index_sims(preds_cod27_sim, area = rep(2*2, nrow(preds_cod27_sim)))
index28_cod_sim <- get_index_sims(preds_cod28_sim, area = rep(2*2, nrow(preds_cod28_sim)))
# Add some columns
index24_cod_sim <- index24_cod_sim %>% mutate(sub_div = "24", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index25_cod_sim <- index25_cod_sim %>% mutate(sub_div = "25", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index26_cod_sim <- index26_cod_sim %>% mutate(sub_div = "26", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index27_cod_sim <- index27_cod_sim %>% mutate(sub_div = "27", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index28_cod_sim <- index28_cod_sim %>% mutate(sub_div = "28", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
### Flounder
# SD 24
preds_fle24_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 211,032 rows (85%), 38,421 rows remaining
# SD 25
preds_fle25_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 178,983 rows (72%), 70,470 rows remaining
# SD 26
preds_fle26_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 184,410 rows (74%), 65,043 rows remaining
# SD 27
preds_fle27_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 224,505 rows (90%), 24,948 rows remaining
# SD 28
preds_fle28_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 198,882 rows (80%), 50,571 rows remaining
# Now calculate the CPUE index (average)
index24_fle_sim <- get_index_sims(preds_fle24_sim, area = rep(2*2, nrow(preds_fle24_sim)))
index25_fle_sim <- get_index_sims(preds_fle25_sim, area = rep(2*2, nrow(preds_fle25_sim)))
index26_fle_sim <- get_index_sims(preds_fle26_sim, area = rep(2*2, nrow(preds_fle26_sim)))
index27_fle_sim <- get_index_sims(preds_fle27_sim, area = rep(2*2, nrow(preds_fle27_sim)))
index28_fle_sim <- get_index_sims(preds_fle28_sim, area = rep(2*2, nrow(preds_fle28_sim)))
# Add some columns
index24_fle_sim <- index24_fle_sim %>% mutate(sub_div = "24", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index25_fle_sim <- index25_fle_sim %>% mutate(sub_div = "25", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index26_fle_sim <- index26_fle_sim %>% mutate(sub_div = "26", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index27_fle_sim <- index27_fle_sim %>% mutate(sub_div = "27", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index28_fle_sim <- index28_fle_sim %>% mutate(sub_div = "28", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
div_index_sim <- bind_rows(index24_cod_sim, index25_cod_sim, index26_cod_sim, index27_cod_sim, index28_cod_sim,
index24_fle_sim, index25_fle_sim, index26_fle_sim, index27_fle_sim, index28_fle_sim) %>%
mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) # convert to tonnes
#> mutate: new variable 'est_t' (double) with 270 unique values and 0% NA
#> new variable 'lwr_t' (double) with 270 unique values and 0% NA
#> new variable 'upr_t' (double) with 270 unique values and 0% NA
Plot the index and save as csv
ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
geom_line() +
facet_wrap(~species, scales = "free") +
geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
geom_line() +
facet_wrap(~species, scales = "free") +
#geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
ggplot(div_index_sim, aes(year, est_t, color = species)) +
geom_line() +
facet_wrap(~sub_div, scales = "free") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
write.csv(div_index_sim, "output/cod_fle_index.csv")
knitr::knit_exit()